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VAR1A.NCE  Ui-.PUCTION  TECHMQUKS  fOR  SUIULATIRC  .‘URKOV  CllAt.VS 


abstract 


Sinulators  frequencly  wish  Co  esclaace  param- 
eters of  Che  limiting  distribution  of  stable  sto- 
chastic processes.  Several  new  methods  for  reduc- 
ing the  variance  of  such  estimates  will  be  proposed 
and  discussed.  The  methods  are  applicable  to  re- 
ge.ierative  Markov  processes  in  both  discrete  and 
continuous  time  as  well  as  to  semi-Markov  process- 
es. The  methods  are  similar  to  the  technique  of 
multiple  control  variables  yet  differ  in  the  im- 
portant respect  that  it  is  not  necessary  to  calcu- 
late the  me.ins  of  Che  controls.  ' This  is  because 
the  controls  are  chosen  in  such  a way  that  their 
means  actually  equal  the  parameter  of  interest. 

The  r-rthods  do  require  a certain  amount  of  com- 
putation to  be  done  before  the  simulation  begins, 
although  their  cost  should  be  relatively  minor  com- 
pared with  that  of  the  simulation.  Numerical  re- 
sults demonstrating  the  effectiveness  of  the  tech- 
niques for  a simple  queueing  model  are  presented. 


I. 


INTRODUCTION 


In  recent  years  computer  simulation  has  become 
a very  important  tool  for  analyzing  the  behavior  of 
stochastic  processes.  As  the  structures  of  widely 
used  processes  become  increasingly  complex  analytic 
results  become  more  difficult  to  obtain  so  that 
frequently  simulation  Is  the  only  computationally 
feasible  method  to  study  a process.  Unfortunately 
simulation  can  be  a very  expensive  tool  to  use.  It 
is  therefore  desirable!  to  develop  methods  that 
can  reduce  the  run  lengths  (and  hence  cost)  of  a 
simulation  yet  still  give  accurate  estimates.  Such 
methods  are  called  variance  reduction  techniq'ies. 
This  paper  will  propose  sevcr.il  new  variance  rediic- 
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tion  techniques  in  the  special  case  when  the  .sto- 
chastic process  being  simulated  is  a Markov  pro- 
cess. 


The  methods  are  all  control  variable  tech- 
niques (see  [5]),  yet  they  differ  from  most  other 
control  methods  in  the  Important  respect  that  the 
means  of  the  controls  need  not  be  explicitly  known. 
This  is  because  the  controls  are  chosen  in  such  a 
way  that  their  means  actually  equal  the  parameter 
of  interest.  The  methods  require  a certain  amount 
of  computation  to  be  done  both  before  and  during 
the  simulation  but  hopefully  their  cost  will  not  be 
so  great  as  to  prohibit  rhe  use  of  the  methods. 

II.  MARKOV  CHAINS 


We  begin  by  Introducing  notation  and  stating 
some  preliminary  results  for  Markov  chains.  The 
reader  should  consult  |1),  12),  or  [3)  fur  a more 
detailed  analysis  of  these  stochastic  processes. 


Let  {X  ,n  > O)  be  a Markov  chain  with  coxint- 
n — 


able  state  space  E = (0,1,2....1  and  n-step 
transition  matrix  p".  We  asstime  tl.e  chain  is 
irreducible,  aperiodic  and  positive  recurrent. 
Under  these  conditions  the  fcllowlng  Proposition 
is  true. 


(1)  PROPOSITION.  There  exlnts  a probability  dis- 
tribution t:  on  E and  a random  variable  X with 
distribution  x such  that 
T,  n 


ij 


Xj  > 0 


for  all  j S E 
• X 


xP.  i.n. 


i 


I 


i-0 


Ij 


(2) 

(3) 

(4> 


The  " " in  (3)  denotes  convergence  in  dislribu- 


tion.  Now  let  fj  be  a real  valued  function  on 


1 


V.irio:uv  KoJuctuin  (confinuod) 


E aaJ  define 


" "f*  ■ I " * r.[f  (X))  (5) 

j J iio  • J J 


It  IS  froqu<*ncly  of  iiitereat  to  know  for  a 


vurivcy  of  functions  f^.  If  the  state  space  is 


very  large  (perhaps  infinite)  tlie  set  of  stationary 
c<;uations  (4)  nay  be  very  difficult  to  solve  so 


that  nust  be  estimated  via  simulation.  It  is 

tlie  efficient  estimation  of  such'  quantities  that  we 


will  concern  ourselves  with.  We  now  develop  an 


alternate  expression  for  r^  which  will  be  useful 


in  simulation. 


Pick  some  sate  in  E (called  the  return 


state)  which  will  be  designated  by  0.  Let  ■ 0 


and  let  T be  the  mth  time  the  process  enters 
m 

state  0 (Tjj  -0).  1. 

j " 0,1,. ..,k  def ine 


state  0 (T_  - 0).  Let  T *7  - T For 

0 m m Q-1 


T -1 
m 


"m(J)  ■ J 

"“'m-1 


(6) 


We  say  that  the  process  is  in  the  nth  cycle  be- 


tween tir.es  T , and  T -1  and  that  T is  the 
m- 1 m n 


length  of  the  mth  cycle.  Because  the  process  re- 


fcaeratos  itself  at  the  times  T (sec  [2))  we  can 

• n 


conclude  that  { Cf  (0),...,Y  (k),  t ).  n > l)  are 
IS  m — 


i.i.d.  (independent  and  identically  distributed) 
random  vectors.  The  importance  of  this  in  a simu- 
lation context  is  that  the  simulation  run  can  then 
be  broken  up  into  randomly  spaced  i.i.d.  blocks  so 
that  the  techniques  of  classical  st  at  i .s  t ics  .can  be 
used  to  analyze  tlie  output  (see  |4)). 


Suppose  now  that  we  simulate  the  process  for 


M cycles.  For  jrach  j define 


M M 

^(M)  - I Y (j)/  I r . 


(7) 


in«l 


m'l 


Uefinc  r and  f(H)  to  lie  the  (k+1)  dimensional 


colu-.n  vectors  witl>  jth  rnlrlcs  r^  and  f ^ (M) 


respectively.  l.et  I bo  the  symmetric  matrix 
whose  (i,f)th  entry  is 


c - EI(Y„(i)  - r,i  )(Y„(j)  - r,:  )1.  (R) 

ij  ffi  inm 


Kc  flssuPC  th.it  c.ich  Plcnjcnt  of  J is  finiti*  md 
tfi)*.  J is  positive  definite,  l.ft  N (0,  ("  ^ ) ) 

dfOAtc-  .1  (^>1)  climensinnnl  r.indo?i  vector  having 


the  niiHivnriate  nornuil  distribution  wi;h  iic.jns  0 


and  J covariance  nintrlx  wiili  cKmcciiCs  O ^ ^ ( t ^ > . 


The  following  pi  opo.si  1 1 on  is  the  key  result  for 
obtaining  point  ostiniufes  and  confidence  interval' 


for  the 


(9)  rROPOSITlON.  If  lT|f  ( < for  all 


J • 0,1,... k and  if  *■“’  for  all  i and 


J ■ 0, 1 , . . . ,k  then 


r_|  - E[V^(J))/E1t^ 


(10) 


fj(M)  "*  r^  a.s.  (almost  surely)  (11) 


>/M(r(M)  - r)  - N (0 , [ /E^ (T ^ ))  . (12) 

As  a corollary  to  (12)  if  we  lot 
e - (e(0) 6(k))  and  let 

0^(3)  ■ 6 )■  e’  (in) 

then  the  following  central  limit  theorem  is  true; 
x/N(ef(M)  - Sr) 


o^(S)/E(;i) 


”*  N(O.l) 


(14) 


where  S(0,1)  denotes  a normally  dS  st  r !hii  i vd 
random  variable  with  moan  0 and  vari.ancc  1. 


111. 


MULTIPLE  F.SriM-rTF.S 


We  now  turn  to  the  first  of  the  variance 
reduction  tecitniques  for  Markov  chains.  Lot  us  fix 
a function  f and  let  r " nf  We  .iro  interested 
in  obtaining  short  confidence  Intervals  for  r. 


To  do  this  we  form  new  .fanctions  f,  so  tiiat 


r^  • rfj  ~ r for  e.ich  J. 
for  f^  let 


As  our  first  candidati 


f , 


f 


for  j - 0,1,. ..,k. 


(IS) 


We  then  have 


Tj  - nfj  - i;(r*f)  - T:r(P'**’f) 

(since  T*  *♦  vP) 


ii(f’"'o 


“ nf 


.5-' 


j-1 


so  that 


r for  all  J.  |iofit\e 


f^(M)  as  in 


(7).  By  tl))  f J (M)  — r for  each  .J.  Now  let 


S he  .some  vector  so  I h.U 


2 


I H(J)  - 1 

J-0 

an.'  1ft  r.(.':)  bn  do(  ined  I'y 

p 

k 

r.(M)  - I B(j)r  CO 
J-0  J 


(16) 


(17) 


Thnn  r .(M)  -*  r a.».  as  M . Now  using  the 

n 

central  lisiii  theorera  in  (14)  aod  oquation  ilt)  we 
have 


\7M(?gC0-r) 

Oj,(B)/e(t^) 


N(0,1). 


(18) 


We  pick  6-8  to  aininite  the  variance  tcro  in 
(IS)  (this  gives  us  Che  shortest  eonddence  Inter- 
val possible).  If  we  let  e be  a (k'tl)  dimen- 
sional row  vector  with  each  entry  equal  to  1 we 
then  h,ive  ; 

6 - e ^ ^/e  I ^e ' 


oJ(3*)  - 1/e  r 


1 . 
c 


(19) 


(20) 


Lot  - 0^(S*)/aQQ.  Combinliifi  f ^(M) , . . . , f|^(H)  in 

the  above  manner  means  that  to  obtain  confidence 
intervals  of  equal  length  we  need  simulate  only 

as  many  cycles  than  we  would  need  if  we  used  no 
variance  reduction  technique. 

This  method  can  be  extended  to  continuous  time 
i;ark.ov  chains  and  semi-Markov  processes  hy  trans- 
forming them  into  discrete  time  MarV.ov  chains  using 
the  techniques  of  (7).  Table  1 lists  the  variance 
reductions  for  cstioaclng  the  expected  queue 
length,  E(X),  in  a finite  capacity  M/M/1  queue. 

We  let  0 - \/\i  where  X is  the  arrival  rate  and 
y is  the  service  rate.  For  case  of  computation 
the  capacity  was  chosen  to  be  14.  It  should  be 
eaphaslred  that  these  figures  are  actual  calcu- 
lations of  variance  reductions  (based  on  the 
methods  of  [7])  and  not  .simtil.it  ion  results. 
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table  1 


ii\nv 

e.  Reductions  for 

Finite  C. 

[M/l 

Oueur  U.sin.((  'Uilc 

i P 1 e E s t i 1 

,.2 

„2 

2 

P 

'^2 

25 

.1163 

.0292 

.0073 

50 

.2341 

.1121 

.0524 

90 

.6050 

.2659 

.1148 

95 

.6880 

. 34  25 

.1607 

99 

. 7404 

.4056 

.204  7 

ip.ic  it  V 


TV.  MATRIX  iterative  TECHNIQUES 

Ue  now  examine  several  other  methods  in  an 
attempt  to  obtain  greater  variance  reductions  for 
a fixed  amount  of  computation  done  before  the  sim- 
ulation. The  basic  idea  here  is  to  partially  solve 
a system  of  linear  equations  (with  some  matrix 
Iterative  technique)  and  then  estimate  the  differ- 
ence between  the  partial  and  true  solutions  via 
simulation. 

We  now  assume  that  E is  finite.  Let  y be 
a vector  with  components  y(i)  defined  by 


y(l)  - E.[  I f(X^)l 
n-0 


(21) 


where  E^I  ) denotes  the  expectation  when 

defined  b 
if  j )•  0 


Xq  “ i.  Let  qP  he  u matrix  defined  by 


"ij 


o^J  - 


^0  if  J - 0 . 

It  can  then  bo  sliown  (see  17))  that 

» 

y-  I or"f 

n-O 

and  that  y s.itisfics 


y - f + qP  y 


(22) 


(23) 


(24) 


Recall  that  r - y(0)/E(t^>.  Let  y"’  be  our 
approximation  to  y after  j iterations  of  the 
Matrix  iterative  method.  We  then  seek  to  find  .s 
function  f^  such  that 


T -1 

y(0)  - y7(0)  + F I \ f^(X  )] 

° n-T  P " 
m-1 


(2b) 


Wc  can  then  set 


T -1 
n 


Vfj)  - y^(0)  + y f fx  ) 
n-T  ^ •' 


(?6> 


TO~1 


3 


anJ  defining  r ^ (M)  as  In  (7)  we  have  r^Cil)  “*  r 
a.s.  Wo  can  tlien  proceed  exactly  nr.  hefore  Co 
obtain  variance  reductions.  Due  to  space  limita- 
tions the  deriv.itions  of  y^  and  for  spocili 

nacrlx  iterative  techniques  are  omitted  here  but 
may  be  found  in  (6).  Table  2 lists  the  calculated 
variance  reductions  for  the  finite  capacity  M/M/1 
using  the  Causs-Seldel  iterative  technique. 


Variance 

TABLE  2 

Reductions  for  Finite  Capacity 

M/M/1 

Queue 

Uslnis  Gauss* 

■Seidel 

„2 

2 

2 

P 

h 

«2 

>^3 

.25 

.0099 

.0015 

.0002 

.50 

.0720 

.0408 

.0207 

.90 

.5939 

.3881 

.1834 

.95 

.6789 

.4  794 

.2414 

.99 

. 7321 

.5446 

.2963 

V.  CONCLUSIONS 

By  viewing  Tables  1 and  2 it  con  be  seen  that 
neither  method  dominates  the  other.  The  variance 
reductions  will,  in  general,  depend  on  the  transi- 
tion matrix  P and  the  function  f.  The  first 
method  has  the  advantage  that  the  variance  re- 
ductions are  independent  of  Che  return  state 
whereas  with  Che  other  methods  care  must  be  taken 
to  pick  a return  state  that  yields  good  variance 
reductions.  It  is  anticipated  that  frequently 
occurring  st.rces  will  produce  the  best  variance 
reductions.  Causs-Seidel  has  the  advantage  that 
any  y®  may  be  used  to  initiate  the  iterative 
procedure.  Ce  erally  speaking  Che  closer 
is  to  y Che  better  the  variance  reduction  will 
be.  This  suggests  that  one  could  simul.ite  a small 
nu.mber  of  cycles  to  obtain  an  initial  y^  and 
then  cor.Tence  the  Gauss-Seidcl  iterations. 

In  order  for  any  of  these  methods  to  be 
computationally  efficient  the  simulator  must  be 
able  CO  compute  (and  store)  the  functions  f ^ . 

The  amount  of  work  involved  in  this  could  be  con- 
siderable unless  the  transition  m.itrix  is  sparse. 

It  is  for  these  types  of  processes  tliac  the  methods 
arc  reco.-.mended . liirther  details  on  the  methods, 
including  more  extensive  numerical  testing,  nay  be 
found  in  16]. 


ACK.hO'..q.KnGKMf::.r . 

The  author  wishes  to  thank  Professors  Donald  L. 
Iglehart  and  Arie  Hotdijk  for  llioir  valu.ible 
suggestions  during  this  research.  The  .luchor  would 
also  like  to  thank  Pat  Rospendowski  for  her  expert 
typing  of  the  nannscripc.  This  research  was 
supported  by  NSF  grant  number  NCS75-23607  and  ONR 
vontract  number  N00014-76-C-0578. 


DIBl.IOCRAPHY 

111  CHUNG,  K.L.  (1967).  Markov  Chains  with 
Stationary  Transition  Probobi 1 1 1 ies,  2nd  edn. 
Springer-Verlag,  Berlin. 

(21  CINLAR,  E.  (1975).  Introduction  to  Stochastic 
Processes.  Prentice-llall , Inc.,  Englewood  Cliff.s, 
N.J. 

(3)  CR^INE,  M.A.  and  IGLEHART,  D.L.  (1974.).  Simu- 
lating stable  stochastic  systems,  II:  Markov 
chains.  J^.  Assoc ■ Comput  ■ Mach.  21,  114-123. 

(41  C?w\NE,  M.A.  and  IGLEHART,  D.L.  (1975).  Simu- 
lating stable  stochastic  systems.  111:  Regenera- 
tive processes  and  discrete-event  simulations. 
Operat ■ Res.  23,  33-45. 

[51  C.4VEK,  D.P.  and  THOMPSON,  C.L.  (1973). 
Programming  and  Probabiliiy  Models  in  Opera! ions 
Research . Brooks/Cole  Pub.  Co.,  Monterey,  Calif. 

|6)  HEIDELBERCER , P.  (1977).  Variance  reduction 
techniques  for  tiie  simulation  of  Markov  processes. 
Ph.D.  Dissertation,  Stanford  University, 

(To  appear) . 

(7)  HORDIJK,  A.,  IGLEHART,  D.L.,  and  SCILASSBERCER, 
R.  (1976).  Discrete  time  methods  for  simulating 
continuous  time  Markov  chains.  Adv.  Appl . Proh ■ 

8,  772-783. 


ii 


UNCLASSIFIED 


SCCuMITv  CLASSIFICATION  OF  THIS  NAOC  >•<»•«  Dim  Sniw*« 


REPORT  DOCUMENTATJON  PAGE 


READ  mSTRUCnONS 
BEFORE  COMPLETING  FORM 


I AI^ONT  numBCN 

41' 


IjToOVT  ACCCUIOH  NO 


S NtClNlCNT-t  catalog  MUMOCN 


A title  rAn<  Suttlilm) 

Variance  Reduction  Techniques  for  Simulating 
Markov  Chains 


S Tvpf  sr  NENONT  • NENIOO  COWENEO 

Technical  Report 


s.  NCMFOmUNO  ONO.  NtNONT  NUMSCN 


I CONTRA^'  oh  OHANT  NUMeCRTO 

N00014-76-C-0578 


7.  Avj^HOAr«> 

Philip  Heldelberger 


10.  NNOONAM  ELEMENT.  NNOJECT  TASK 

ANEA  A NONK  unit  nUMSCRS 

(NR  042-343) 


s nenfonmino  onoanization  name  ANO  AOODESS 
Department  of  Operations  Research  ^ 
Stanford  University 
Stanford,  CA  94305 


II.  CONTNOLLINQ  OFFICE  NAME  ANO  AOOItESS 

Statistics  and  Probability  Program 
Office  of  Naval  Research  (Code  436) 
Arlington.  Virginia  20360 


IS.  NE^NT  DATE 

September  1977 


IS.  numser  OF  FAoes 


14.  MONITONINO  AGENCY  NAME  0 A00NBSS</(  dlllmtmn  /fF*  Ca0Ut1llng  OIHf) 


IS.  security  class.  (•!  lAla  raFWt) 

UNCLASSIFIED 


ISa  OECLASitFICATION/DORHORADlNO 
' ISDULE 


SCnI 


IS  OISTNIEUTION  statement  (oI  IMt  R«F«r«) 

Approved  for  public  release:  distribution  unlimited. 


IT.  OISTRISUTION  STATEMENT  ,'a<  CDm  tbtirtt  MttMFS  In  Block  SO,  II  ^llottnl  Irom  MoBOtt) 


At  th 


IS./SURRLEMEMTARV  notes 

This  paper  will  be  presented  At  the  1977  Winter  Simulation  Conference, 
December  5-7,  1977,  National  Bureau  of  Standards,  Caithersburg,  Maryland. 


It  <CV  VOffOl  (Comtimaf  it  w— m4  hr  himcM  mmhrt) 


SIMULATION,  VARIANCE  REDUCTION,  MARKOV  PROCESS 


aO.  A9STI9ACT  (Cmthmm  «m  «!#•  ti  hr  hl—h 


see  attached 


UlSClASS  IFIED 


ttCuMiTv  Ck«tti'iCATiON  or  This  r«ac  Oms  MttttMO 


Variance  Reduction  Techniques  for  Simulating  Markov  Chains, 
by  Philip  Heidelberger 

Simulators  frequently  wish  to  estimate  parameters  of  the  limiting 
distribution  of  stable  stochastic  processes.  Several  new  methods  for 
reducing  the  variance  of  such  estimates  will  be  proposed  and  discussed. 

The  methods  are  applicable  to  regenerative  Markov  processes  in  both  discrete 
and  continuous  time  as  well  as  to  semi-Markov  processes.  The  methods  are 
similar  to  the  technique  of  multiple  control  variables  yet  differ  in  the 
important  respect  that  it  is  not  necessary  to  calculate  the  means  of  the 
controls.  This  is  because  the  controls  are  chosen  in  such  a way  that  their  | 
means  actually  equal  the  parameter  of  interest . The  methods  do  require  a 

I 

certain  amount  of  computation  to  be  done  before  the  simulation  begins,  j 

although  their  cost  should  be  relatively  minor  compared  with  that  of  the 
simulation.  Numerical  results  demonstrating  the  effectiveness  of  the 


techniques  for  a simple  queueing  model  are  presented. 


ItCuetTV  CkMWriCATION  Of  TMII  f ASinnkMi  Oats 


